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Abstract 


A finite element code is developed for solving 3-D, dynamic, elasto-plastic, contact 
problem, including the effects of friction at the contact interface. Elasto-plastic behaviour 
is modelled by an associated flow rule based on the von Mises yield criterion and power law 
type isotropic hardening. Jaumann stress rate tensor and incremental Green-Lagrange 
strain tensor are used as the objective stress and non-linear strain measures respectively. 
The material and geometric non-linearities are handled by using the updated Lagrangian 
method. A node-to-segment interface model is used for developing the contact stiffness 
matrix. Contact constraints are imposed using Lagrange multiplier method. Coulomb’s 
friction law is used for calculating the friction forces. A finite element - finite difference 
scheme is used for spatial and temporal discretization respectively and Modified Newton- 
Raphson iteration method is used for solving the non-linear governing equation. To 
improve the convergence characteristics of Newton-Raphson iterations, under-relaxation 
technique and line search method are used. The finite difference scheme used for time 
integration is Newmark’s algorithm. Unloading is incorporated to include the effects of 
global or local unloading during the analysis. Skyline scheme of assembly and static 
condensation scheme are used to reduce the computational resources required for the 
analysis. The finite element code is validated using three standard problems from the 
published literature. Finally, a 3-D, dynamic, elasto-plastic contact problem with friction 
is solved to demonstrate the applicability of the code. 
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Chapter 1 


Introduction 


Contact interactions occur in virtually all structural and mechanical systems. They can 
be intentional as in the case of transfer of loads between two components or unintentional 
as in a vehicle crash. If the contacting bodies are not in static equilibrium during the 
contact, then it is called a dynamic contact or impact. The basic characteristics of a 
dynamic contact phenomenon are very brief duration, high force levels reached, rapid 
dissipation of energy and large accelerations and decelerations present. By nature, the 
contact phenomenon always involve friction. In many cases, the bodies undergo plastic 
deformations during contact. A rigorous analysis of the contact problem requires the 
consideration of all these complex aspects. The analysis of the problems like the autmobile 
crashes involves the use of dynamics, elasto-plasticity and large deformation in a contact 
analysis scenario. 

Closed form mathematical solutions are extremely difficult to develop even for simple 
contact problems. As the contact scenario becomes more complex, as in the case of elasto- 
plastic materials, arriving at a closed form solution is nearly impossible. This necessitates 
the use of approximate numerical tools for the analysis, of which finite element analysis 
is arguably the most versatile option. . 

1.1 Review of Literature 

The classical incremental theory of elasto-plasticity [1] is generally taken to be the basis 
for predicting elasto-plastic response. The constitutive equation for these materials are 
derived by combining the stress-strain relations for elastic and plastic response. The 
constitutive relation for a plastic response requires three fundamental characteristics of 



the material: a yield criterion, a flow rule and a strain hardening relationship. The yield 
criterion demarcates the elastic behaviour from the plastic behaviour. It specifies the 
state of multi-cixial stress corresponding to the start of the plastic flow. Various yield 
criteria have been proposed by Tresca, von Mises etc. [2], where von Mises criterion is 
found to be the simplest and the most accurate in terms of the predicted values. The 
flow rule relates the plastic strain rate to the current stress. It is represented in terms 
of a plastic potential, analogus to the elastic potential, and is termed as the associated 
flow rule if the potential is taken as the yield function. The hardening rule specifies how 
the yield condition is modified during the plastic flow. For the purpose of analysis, the 
non-linear hardening relationship is usually modelled in terms of a power law. 

Unlike the case of small deformation of elastic materials, the governing equations 
of elasto-plastic materials, which undergo large deformation, contain two t5rpes of non- 
linearities: (i) geometric non-linearity .associated with non-linear strain displacement re- 
lations and (ii) material non-linearity in the form of non-linear stress-strain relations. An 
incremental formulation like the Updated Lagrangian Formulation is the most convenient 
way to handle such non-linearities unlike the total Lagrangian formulation. The advan- 
tage of the former is that, if the increment size is sufl&ciently small, linearized incremental 
equations can be used. However, when the overall deformation is large, the number of 
increments is to be prohibitively large. In such a situation, we have to resort to some 
iterative scheme for solving the non-linear incremental equations. The prominent among 
such schemes is the Newton-Raphson method. In this method, the stiffness matrix is up- 
dated in each iteration of the increment. Thus, we work with the tangent modulus at 
that instant for each iteration. This method has two shortcomings. The first being that 
there is chance of divergence if the slope of force-displacement curve increases with the 
increment. Secondly, it takes more computational time as the stiffness matrix has to be 
computed afresh in each iteration [3], Thus, the Modifled Newton-Raphson method in 
which, the stiffness matrix is kept constant during the increment, is adopted. 

The dynamic analysis is very basic to studying impact problems. There are various 
ways of performing the time discretization in this category of problems of which, the 
direct integration methods are prominent. These methods include the central difference 
method, Houbolt method, Wilson-0 method and Newmark method [4]. The first method 
is conditionally stable one, while the rest of them are unconditionally stable. Analysis 
of accuracy and speed of different algorithms is presented by Bathe [4]. It has been 



concluded that although central difference method is beneficial in terms of speed and 
cost, the trapezoidal method of Newmark’s algorithm is best in terms of stability and 
accuracy. 

1.1.1 Elasto- plastic Finite Element Analysis 

As stated earlier, the classical theory of elasto-plasticity [1] is generally taken to be the 
basis for predicting the elasto-plastic response. Argyris [5] presented one of the earliest 
formulations and solution techniques for elasto-plastic problems. The effect of plastic- 
ity was simulated by “initial loads” and the matrix method was used to discretize the 
linearized governing equations. Two approaches, the direct incremental method and the 
iterative incremental method, were discussed for obtaining the solution. The eaxly appli- 
cations of the iterative incremental method [3] used the initial stiffness method. However, 
this approach is not suited when the overall deformation is very large as it results in very 
slow convergence, or in many cases, even divergence. 

Nagtegaal and Jong [6] dealt with the computational aspects of elastic-plastic large 
deformation analysis. A modified set of governing equations were formulated, which were 
applicable for small strain and large rotational increments. Special attention was paid to 
the integration of these equations in the deformation history. It is predicted that explicit 
methods, like the tangent modulus procedure, are conditionally stable whereas the implicit 
schemes, such as the mean normal method, are preferable. The main difficulty with the 
three dimensional analysis is the vast computational memory that it requires. But, with 
the advent of better machines and evolution of more efficient algorithms this problem has 
been solved and more and more such analyses are being acrried out. 

Since the elasto-plastic constitutive relation is generally given in the differential form, 
the constitutive function must be integrated to obtain the incremental stress. The pro- 
cedures commonly used are the simple and robust Euler forward integration technique, 
the implicit Euler backward integration technique and the radial return method. The 
trend has lately been shifting towards the use of consistent schemes which improve the 
convergence characteristics of the equilibrium iterations [7]. 

Bathe [4] has briefly reviewed the different approaches used for obtaining the solution 
of general non-linear problems. Total and updated Lagrangian descriptions of motion 
have been used to formulate the incremental equations of motion. It has been noticed, for 



the problems solved, that the two formulations give identical results, provided appropriate 
constitutive relations are used. One of the areas of controversy in analyzing large defor- 
mation is the choice of appropriate objective stress measures. Various objective stress 
measures have been discussed by Bathe [4] and Chakrabarty [8]. The classical theory of 
elasto-plasticity assumes that the incremental strain tensor can be additively decomposed 
into the elastic and plastic parts. While this may be true for small incremental deforma- 
tion, this assumption is no longer valid for large incremental strains and rotations, unless 
appropriate strain measures are used [9, 10] . 

1.1.2 Dynamic Finite Element Analysis 

Bathe et al [11] were probably the first to propose a large deformation dynamic finite 
element formulation. Total and updated Lagrangian descriptions of motion were used 
to formulate the incremental equations of motion. A detailed finite element analysis of 
non-linear static and dynamic response was also carried out by Mondkar and Powell [12]. 
The purpose of this paper was to review the theoretical and computational techniques 
for non-linear structural analysis, and their applications in a general computer program. 
A Lagrangian frame of reference was used and a variational form of the incremental 
equations of motion for large displacement, finite strain deformation was incorporated in 
the finite element program. It was found that no specific time integration scheme was 
optimal for all types of non-linear behaviour. It was inferred that for large load steps it 
was essential to use iterations to obtain accurate results. 

Direct integration methods are the most commonly used methods for solving the dis- 
cretized and linearized governing equations of motion. Dynamic equilibrium including the 
effect of damping forces is sought at discrete time instants instead of at all times. The 
advantage is obvious: all solution techniques employed in static analysis can probably be 
used most effectively in direct integration. The selection of any one direct integration 
methods depends upon consideartions of accuracy, stability and cost. The importance of 
using equilibrium iterations in dynamic non-linear analysis was first realized by Bathe [4]. 
An extensive discussion of these and other methods can be found in Bathe [4]. 

More recently, Gendy and Saleeb [13] found that, for a more comprehensive perspec- 
tive, work on the non-linear dynamic responses of plates and shells calls for detailed 
studies of several important factors. These include the effect of large spatial rotations on 



the geometric stiffness and inertia operators, accurate updating procedures for nodal rota- 
tions and associated angular velocities and accelerations, as well as material inelasticities 
(especially for finite strains). Several of these issues were examined in conjunction with 
a newly developed mixed finite element formulation for plates and shells. To this end, 
and restricting the scope to the case of large overall motions but small strains, low-order 
displacement/strain interpolations are utilized, together with a radial return algorithm 
(backward Euler-integration scheme) for plasticity effects. The Newmark implicit scheme 
has been employed to intgrate the semi-discrete equations of motion. 

1.1.3 Contact Analysis 

The analysis of contact has been a fascinating area since the initial days of the study of 
mechanics itself. Many investigators tried to determine the stresses and displacements in 
two elastic bodies coming in contact. In three dimensions, the problem of two contact- 
ing elastic bodies was first formulated and solved by Hertz [14] under several restrictive 
assumptions. Hertz assumed that the geometry of general curved surface in the vicinity 
of contact can be described by quadratic terms only. Neglecting friction and assuming 
that the bodies in contact deform as though they were elastic half-spaces. Hertz used 
the Boussinesq solution to deduce that the contact pressure distribution has to be el- 
liptical to produce an elliptical area. For the contact between two spheres, the contact 
area is circular. Research work based on Hertz’s theory up to 1982 has been reviewed 
by Johnson [15]. Although Hertz evaluated the dimensions of contact area and stresses 
at the contact surface, he presented only a speculative sketch of the subsurface stresses. 
The stresses beneath the contact surface were analyzed later. For circular contact area of 
radius a, it was observed that the maximum shear stress occurs at a depth of 0.57a [15]. 
Later on Muskhelishvili [16] and Gladwell [17] solved some small deformation elastic con- 
tact problems using the analytical techniques. 

With the development of computers, people resorted to numerical methods. Finite 
element method has been the most widely used. In earlier studies, contacting bodies 
were assumed to consist of linearly elastic materials and hence were assumed to undergo 
only small deformation. The loading was assumed to be proportional. The node-to- 
node contact interface model was used for discretizing the contacting bodies and the 
contact conditions (constraints) involving nodal displacements and forces were applied 



by modifying the combined stiffness matrix or the flexibility matrix of the contacting 
bodies. The investigators who used this model are Ohte [18], Gaertner [19], Prancavilla 
and Zienkiewicz [20], Sachdeva and Ramakrishnan [21] etc. 

Kalker and Randen [22], Hung and de Saxce [23] and Mahmoud et al. [24] used the 
technique of mathematical programmming for solving the frictionless small deformation 
elastic contact problem. They minimised the total strain energy to find the contact area 
and the contact pressure distributions. 

When the deformation is large, there is possibility that a pair of nodes presently in 
contact, may undergo a large relative displacement. If a node-to-node interface model is 
used, re-meshing of the bodies becomes necessary for carrying out the subsequent analysis. 
To avoid frequent re-meshing, one must use the node-to-segment interface model, in which 
one node of a contacting body comes in contact with a segment of the other. When a 
node-to-segment interface model is used, the contact conditions acquire a complex form, 
when expressed in terms of nodal variables. To impose these conditions, an additional set 
of finite element equations, involving the contact stiffness matrix, needs to be developed 
using the principle of virtual work of the contact forces. The methods, which are normally 
used for implementing the contact constraints are Lagrange multiplier method, penalty 
funtion method and transformation matrix method. Bohm [25] presented a comparison of 
different contact algorithms with applications. A survey of contact algorithms is given by 
Bourago [26]. 

Kikuchi and Oden [27] have provided the mathematical foundation for the Lagrange 
multiplier method. The work of Bathe and Chaudhary [28] seems to be the first at- 
tempt to develop an expression for contact stiffness matrix using this method. This was 
done for planar and axisymmetric problems assuming the contact segment to be straight. 
Lagrange multiplier was used to impose the contact constraints including friction forces 
based on Coulomb’s friction law. Non-linear contact kinematics was used for developing 
a consistant contact matrix for two dimensional problems by Wriggers and Simo [29] and 
by Parisch [30]. They used both the Lagrange multiplier method and penalty function 
method. The transformation matrix method for two dimesional problems was proposed 
by Chen and Yeh [31]. A formulation for non-linear frictional model was proposed by 
Gallego and Anza [32] using perturbed' Lagrangian functional. All the above formulations 
were used for static large deformation elastic contact problems. Some of these formula- 
tions were extended to elasto-dynamic contact problems [4, 29] and dynamic elasto-plastic 



problems [33]. A literature survey of contact dynamics modelling is presented by Gilardi 
and Sharf [34] . Sung and Kwak [35] and Bittencourt and Creus [36] present the analysis 
of dynamic, large deformation contact including the effects of friction. Master’s thesis by 
Werner [37] is an attempt to compare 'the experimental results and the results predicted 
by a commercial finite element software for an elasto-plaatic impact problem. 

Contact-impact algorithms, which are sometimes called as slideline algorithms, are the 
computationally time-consuming parts of many explicit simulations of non-linear prob- 
lems, becuase they involve so many branches, and so are not amenable to vectorization, 
which is essential for speed on supercomputers. Belytschko and Neal [38] introduced the 
pinball algorithm, which is a simplified slideline algorithm and is readily vectorized. Its 
major idea is to embed pinballs in surface elements and to enforce the impenetrability 
condition only to these pinballs. It can be implemented either by a Lagrange multiplier 
or penalty method. Malone and Johnson [39] developed a parallel contact algorithm and 
implemented it in a non-linear dynaniic explicit finite element program to analyze the 
three-dimensional shell structures. The contact algorithm accounted for initial contact, 
sliding and release through the use of parametric representation of motion of points lo- 
cated on the surface of structure combined with a contact surface representation, which 
approximates the actual surface by means of triangular search planes. Brown et at [40] de- 
scribed a general strategy for parallelizing solid mechanics simulations. Such simulations 
often have several computationally intensive parts, including finite element integration, 
detection of material contacts and particle interactions if smoothed particle hydrodynam- 
ics is used to model highly deforming materials. Their strategy is to load-balance each 
of the significant computations independently with whatever balancing technique is most 
appropriate. The main advantage is that each computation can be scalably parallelized. 
The drawback is the data exchange between the processors and the extra coding that 
must be written to maintain multiple decompositions in a code. 

The problems of undesirable oscillatory solutions and stability of time-integration 
schemes have been tackled by quite a few people. Taylor and Papadopoulos [41] address 
the problem of oscillatory solution along the contact interface. The standard second order 
implicit integrators of the Newmark family, used in earlier works, were examined by them. 
In their work, control of these oscillations is attempted by introducing an artificial bulk 
viscosity for the compressive waves in each body. Laursen and Chawla [42] proposed an 
energy conserving algorithm to control the oscillations for the frictionless dynamic contact 



problems. It is shown that a Lagrange multiplier enforcement of an appropriate contact 
rate constraint produces these conservation properties. In particular, the ability of the 
formulation to produce accurate results, where more conventional integration schemes 
fail, is emphasized by numerical simulations. Armero and Petocz [43] presented a new 
dissipative time-stepping algorithm to improve the stablity of time integration scheme. 
They showed that although the traditional time-integration schemes were stable for most 
of the linear problems, instability would creep in as soon as we try to use them in non- 
linear cases. Given these consideration^!, they developed implicit time-stepping algorithms 
for contact problems that posess unconditional (energy) stability in time and lead to stable 
enforcement of contact constraints. 

Quite a few researchers have proposed new techniques. Ayari and Saouma [44] devel- 
oped a new formulation for contact problems using the concept of fictitious forces. Con- 
trary to the most existing models, this one involves very few matrix decompositions and 
thus is computationally inexpensive. The model was applied to fracture mechanics based 
analysis of a cracked dam under seismic excitation. They also showed that their model 
could be applied to a vast majority of contact problems and was not only restricted to 
fracture phenomenon. Simo and Laursen [45] extended the augmented Lagrangian fomula- 
tion, which has certain advantages over traditional Lagrangian formulation, to problems 
involving frictional contact. Kane et al [46] proposed a nonsmooth contact algorithm, 
which is capable of dealing with multibody nonsmooth contact geometries for which nei- 
ther the normals nor the gap functions can be defined. Farahani et al [47] presented a 
solution method for contact/impact problems, which is different from Lagrange multiplier 
and penalty methods. This method is based on the stiffness matrix transformation and 
eliminating the normal degree of freedom of contactor node. The method is absolutely 
general and can be used in static as well as dynamic non-linear problems. In this case, 
the number of unknowns does not increase and the solution scheme needs no modification 
as is the case with prior methods. Besides, the boundary conditions are satisfied exactly. 

1-2 Objective of the Present Work 

■lir}i e obf ective of the present work is to develop an efficient finite element software for 
d(’<formation, dynamic, elasto-plastic contact problem including the effects of 
(^"Te contact interface. Since we have large deformations, including slipping 



at the contact interface, a node-to-segraent contact model is used in contact analysis. 
A method outlined by Zhong [33] is used for calculating the contact stiffness matrix. 
Lagrange multiplier method is used for applying the contact constraints as penalty method 
is found to be dependent on the penalty number, which is a free (user defined) parameter. 
Coulomb’s friction law is used for modelling the frictional effects at the contact interface. 

Elastic behaviour is modelled using generalised Hooke’s law and the plastic behaviour 
is modelled by an associated flow rule based on von Mises yield criterion and power- 
law type isotropic hardening. The total strain is obtained by adding the elastic and 
plastic components. Jaumann stress rate tensor and incremental Green-Lagrange strain 
tensor are used as objective stress and non-linear strain measures in the formulation of 
constitutive equations. Incremental analysis is carried out using the updated Lagrangian 
method [4]. Modified Newton- Raphson iterative method is used for solving the resulting 
finite element equation. In each Newton-Raphson iteration, contact iterations are carried 
out to find the contact reactions. A finite difference scheme is used for carrying out the 
discretization in time. Newmark’s algorithm, which is found to be the best in terms 
of stability and accuracy among the various finite difference schemes, is used for this 
purpose. Unloading scheme is incorporated to take care of the local unloading during 
the analysis. To reduce the computational resources required for performing the analysis, 
static condensation of stiffness matrix for contact and non-contact nodes and skyline 
scheme of assembly are used. Divergence handling techniques, viz., under relaxation and 
line search, are used for improving the convergence characteristics of Newton-Raphson 
iterations. 


1.3 Structure of the Thesis 

Chapter 2 deals with the mathematical formulation of the governing equations for dy- 
namic, large deformation, elasto-plastic problems. The development of the finite element- 
finite difference scheme based on these governing equations and the numerical scheme 
used are discussed in chapter 3. Chapter 4 presents the contact formulation, wherein the 
formulation of the contact stiffness matrix and the calculation of contact reactions are dis- 
cussed. Chapter 5 is devoted for the results of the validation problems used for validating 
the code and an additional problem to show the capabilities of the finite element code. 
The last chapter summarises the work carried out and the scope for future developments. 



Chapter 2 


Mathematical Modeling of Dynamic 
Blast o- Plastic Problem 


In a dynamic elasto-plastic problem, the deformations and strains in the body no longer 
remain small and hence the response of the body is not proportional to the load. The 
response of the body becomes non-linear. The non-linearities can be broadly classified 
into two: Geometric non-linearity and material non-linearity. Both these non-linearities 
occur in the analysis of large deformation of elasto-plastic materials. The most convenient 
formulation for such analysis is the Updated Lagrangian Formulation with the constitutive 
relations represented in an incremental manner. 

2.1 Introduction 

In the study of the deformation of a body subjected to external loading, often the original, 
undeformed and unstressed state of the body is used for the formulation of its governing 
equation. This method is known as the Lagrangian formulation. This formulation is 
convenient for small deformation problems, where the deformed configuration does not 
deviate from the original configuration significantly. Hence, the deformations are pro- 
portional to the loading and can be described by an infinitesimal or linear strain tensor, 
for which the strain - displacement relations are linear. On the other hand, for large 
deformation problems, one has to use a finite strain measure, which is expressed by a 
non-linear strain displacement relation. Furthermore, the equations of motion, expressed 
in the original configuration, depend on the deformation. This makes the governing equa- 
tions highly complex and too difficult to solve. In such cases, it is convenient to solve 



the problem in an incremental manner known as the Updated Lagrangian Formulation. 
In this formulation, the state of the body is updated in increments using the calculated 
incremental deformations and stresses. The response during At is calculated based on 
the configuration at time t, where At represents the increment in time, called time step. 
Hence, all the relations are to be in an incremental form rather than the total form. This 
methodology is particularly useful for elasto-plastic materials, since the constitutive rela- 
tions for these materials depend on the state of the body and therefore they are expressed 
in an incremental form. 


2.2 Stress Measures 


The large deformation problems are normally associated with rigid body rotations. Hence 
the stress measures used in these problems should be objective, which means that they 
should vanish if the increment is a pure rigid body rotation. The most commonly used 
stress measure known as Cauchy stress tensor is not objective and hence it can not be 
used in a constitutive equation. There are numerous objective stress measures, of which 
two are discussed below. 

One of the commonly used objective stress measures is the second Piola-Kirchoff stress 
tensor which is related to Cauchy stress tensor using the concept of 

equivalent work between the two configurations [4]: 


t+AtQ 




t t+At 

t+At ^mn 




( 2 . 1 ) 


Here */? and represent the densities at time t and time t + At respectively, t+At 
is the derivative of the position vector at time t, viz., ^Xi with respect to Another 

frequently used objective stress measure is Jaumann’s stress rate tensor aij. It is related 
to the Cauchy rate dy by the relation [4]; 


‘(7y At = &ij At - hikCQjk At) - ^ajkCQ,ik At) (2.2) 

fijfc At = — (tAujj tAuj^i) (2.3) 

Here, ‘Oy* At represents the incremental spin tensor and Auij denotes the derivative of 
the incremental displacement vector tAuj with respect to hsj. 



The second Piola-Kirchoff stress tensor is energy conjugate to the Green-Lagrange 
strain tensor. Therefore, it is used in this work to develop the virtual work expressions in 
Updated Lagrangian formulation, whereas for developing the elasto-plastic constitutive 
equations, the Jaumann stress rate tensor, which is closely related to Cauchy stress rate, 
is used. 


2.3 Strain Measures 


The linear or infinitesimal strain tensor, normally used in the small deformation problems, 
is not a valid choice for the large deformation, elasto-plastic analysis. Green-Lagrange 
strain tensor is one of the most commonly used strain measure for large deformation 
analysis. 

The incremental Green-Lagrange strain tensor is a non-linear function of the displace- 
ment [4]: 

(Aejj = — ((Aujj H- fAzij^j -f- tAu^jj) (2'4:) 

When the incremental deformation is small, the incremental linear strain tensor 


jASy 


-(fAuj-j -|- fAuj^j) 


(2.5) 


can be used as the measure of deformation. 


2.4 Elastic Constitutive Equation 

The incremental elastic constitutive equation is the generalized Hooke’s Law relating the 
increment in Jaumann stress components and the elastic part (tAe? ) of the incremental 
Green-Lagrange strain components: 

iAcTy = C^-^i (Ae^j (2.6) 

The tensor for isotropic case is given by: 

= 2ia5ik5ji (2.7) 


where A and // are Lame’s constants. 



2.5 Elasto-plastic Constitutive Equation 


As the stresses in the material exceed the yield stress, the material undergoes plastic 
deformation. The elastic constitutive relation described above is not applicable for the 
plastic deformation. An elasto-plastic relationship between the incremental stress and 
strain tensors based on the von Mises yield criterion and isotropic hardening is developed 
below. 

For an isotropically hardening material, the plastic potential is given by [1]^: 

F'{<^ij,p) =o-eg{(rij) - CXyip) ( 2 . 8 ) 


where, 


F = 0 


(2.9) 


represents the yield criterion. The plastic potential F depends on the Cauchy stress tensor 
through the second invariant of its deviatoric part, called as equivalent stress, and defined 
by: 


^eq — 



( 2 . 10 ) 


where ah is the deviatoric part of Uij. Further, F depends on the variable yield stress of 
the material, ay, through the hardening parameter p. Assuming that the strain hardening 
hypothesis is valid, p is the equivalent plastic strain . Therefore, 


V = ely = j(klq ( 2 . 11 ) 

where, 

(I <&»,■*«)" (2-12) 

Here, deh is the plastic part of the incremental linear strain tensor and the integration 
in equation (2.11) is to be carried along the paxticle path. The dependence of ay on p is 
approximated by a power law type relationship: 


a, - (a,)o = K{el,r (2.13) 

Here, (aj,)o is the yield stress at zero plastic strain, K is the hardening coefficient and n 
is called as the hardening exponent. 


Hn this section, the superscript/subscript denoting time have been omitted for convenience 



The plastic part of incremental linear strain tensor de^j is obtained from the plastic 
potential using the following relation: 

dF 


d4j = dX 
^ oau 


(2.14) 


y 


where dX is a scalar. This equation is called as Flow Rule. Differentiation of equation (2.8) 
with respect to <7^ gives: 


dF 

dvij 




Therefore, dX is determined as: 


d\ = ckf. 


(2.15) 


(2.16) 


Further, the hardening relationship and the yield condition can be used to express dX as: 


dX - 
dX- — 


da, 


eq 


H 


where, 


= ^ = KnK)"-' 


(2.17) 


(2.18) 


is the slope of the hardening curve. Substitution of equations (2.15) and (2.17) in equa- 
tion (2.14) leads to the following constitutive equation: 


de^ = ^ rr' 

^ 2Haeg 


(2.19) 


This constitutive relationship between the deviatoric stress tensor and the plastic part 
of the incremental linear strain tensor is not really convenient for the Updated Lagrangian 
formulation for which an incremental stress-strain relationship is needed. This can be 
obtained from equation (2.19) as follows: 


Serb daeg. 

Note that, from equations (2.8) and (2.15), we get: 


da. 


eg 


dF 


daki daki 2a, 


'^kl 


(2.20) 


(2.21) 


eg 


Substitution of equation (2.21) in equation (2.20) leads to the following incremental plastic 
stress-strain relationship: 


d£% 


ATT 2 

^Halg 


(2.22) 


The incremental elastic stress-strain relationship (equations (2.6) & (2.7)) in terms of 
the incremental linear strain can be expressed as: 

detj = ^[-J^dakkdij + (1 + i^) day] 


(2.23) 



where, is the elastic part of E is the Young’s modulus and u is the Poisson’s ratio. 
Adding the two relationships, we get: 


£ij 


p®. -J_ p?. — 

^%3 ^ ^13 


-g-%fe+ -^S,t6j,+ 


eg J 


da. 


kl 


(2.24) 


This is the incremental elasto-plastic stress-strain relationship needed in the updated 
Lagrangian formulation. However, it is the following inverse relationship which is more 
useful: 


da, 


ij 


^fjkl ^kl 


(2.25) 


where, 

Here, ji is one of the Lame’s constants (also called as shear modulus) and v is Poisson’s 
ratio. 

Note that the stress increment appearing in equation (2.26) must be an objective stress 
increment in the sense that daij must be a zero tensor if the increment is a pure rotation. 
The incremental stress measure used in the present work is the Jaumann stress measure 
discussed in section 2.2. 

The relationship (2.25) has been derived assuming the increment size to be small 
and using the incremental linear strain tensor as the strain measure. When a large 
size increment is to be used along with the incremental Green- Lagrange strain measure 
(defined in section (2.3)), the derivation is similar. The relationship (2.25), when the 
stress and strain measures are replaced respectively by the incremental objective stress 
measure and the incremental Green-Lagrange strain measure, takes the following form: 

. ^t+At 

= l ‘C^SdiAem) (2.27) 


where. 


^C^jki = 2 ^^ Sik Sji + 


p 


-did Si 




1- 2(3)U+ Kn{^e’^q)^~^)*ajgJ (2-28) 

Here, H has been replaced by the expression (2.18) and the left subscript t has been added 
to make it explicit that these quantities are to be evaluated at time t. 


2.6 Incremental Updated Lagrangian Formulation 

The objective of the updated Lagrangian formulation, is to establish dynamic equilibrium 
in the configuration at time t-\- At when all the static and kinematic variables at time t 



are known. The principle of virtual work requires that [4]: 


/. 






% 5C+^%) d^+^V + (2.29) 


where, 

= density at time t + At 
— acceleration at time t + At 

_ virtual displacement at time t + At 
t+Aiy _ volume at time t + At 

= Cauchy stress tensor at t + At. 

The virtual linear strain is defined as: 


Further, 



(d5*+^% 


d^+^% ) 


(2.30) 


t+AtR 


L 
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jt-^Aty 

f 

t+^iSc 


i+A^. ^(t+A<^.) ^+Ay^ 


I 


t+Ai 5 „ 


(2.31) 


where, 

= body force per unit volume at time t + At 
‘+'^^ti )5 = applied traction per unit area at time t + At 
= surface at time t + At with traction specified 
^'^^X'^i)c — contact traction per unit area at time t + At 
= contact surface at time t + At 

The main difficulty in applying equation (2.29) is that the configuration at time t+ At 
is unknown. Therefore, the virtual work expression at time t + At is transformed to 
an integral over the volume at tim t by using the principle conservation of mass. It is 
assumed that the external load term (2.31) is deformation independent for the formulation 
of governing equation. The expression (2.29) after the transformation becomes [4]: 


f V d^V+ f d^V= ‘+^‘1? 

Jty JiV 


where the virtual green-Lagrange strain tensor is defined by: 



'dS*+^hii 


8*+^% ^t+A^^, j 


(2.32) 


(2.33) 



The second Piola-Kirchoff stress- tensor can be decomposed as; 


tSij — ISij + t^Sij — ^cTij + jAiSy 

Since, 

the virtual Green- Lagrange strain tensor can be decomposed as: 

= (^GAcy) = S(tAeij + tA^ij) 

where, 

tAsij — GAujj -|- fAuj^i) 

tATjij — ((A'Ufcy- ^Auj-j') 

Therefore, the equation (2.32) can be written with incremental decomposition as: 

V 6(tAui) d^V+ tASij SitAsij) d^V + 

[ tASi, SitAri,j) d^V + f 6{tAvij) d^V + 
jty J^y 

Vij- SitAsij) d*V = (2.39) 

The above equation is simplified by neglecting the third integral, which is a higher order 
term, and approximating tASij as tAski- 

V 5{tAui) d^V+ tAski SitAsij) d^V + 

% 6{tAr}ij) d^V + *aij SitAsij) d^V = (2.40) 

The linearized equation, when solved, will yield only approximate displacement, ve- 
locity, acceleration, strain and stress fields. The approximate quantities are denoted by 
a right superscript (1). The error due to the approximation involved is calculated from 
equation (2.29) as: 

Error = f - 

f t+AUl) t+At^(l) ^t+Ayil) (2.41) 

Jt+AiyO) 

This error is generally minimised by a predictor - corrector scheme. 


(2.34) 

(2.35) 

(2.36) 

(2.37) 

(2.38) 



Chapter 3 


Finite Element-Finite Difference 
Formulation 


The mathematical formulation of the Dynamic Elasto-plastic problem is given in the equa- 
tion (2.40) of chapter (2). However, this equation contains volume and surface integrals, 
which can not be solved analytically for most of the physical problems. In many of the 
cases of interest, the volume and surface may not be expressed mathematically. This ne- 
cessitates the use of some numerical methods for the solution of the governing equation. 
Since it approximates the field variable and geometry simultaneously, the Finite Element 
Method is probably the most powerful method for solving the complex boundary value 
problems. The integration in the time domain, which is an initial value problem, is carried 
out efficiently using a Finite Difference Scheme. This chapter presents the development 
of a finite element - finite difference scheme for solving the linearized governing equation 
of the dynamic elasto-plastic problem. A numerical iterative technique for minimizing the 
error term given by equation (2.41) is also presented. 

3.1 Matrix Notation 

Matrix notation is used in the development of finite element - finite difference equations 
for the convenience in computer implementation. 

The components of the strain tensors jAsy and are represented in the array 
form as follows: 


({As}^ — {tASxxj tASyy, ^ ^tAs^z} 


(3.1) 



^|A77} • • •} 


(3.2) 


Note that the shear strain components in the array i{Ae}^ contain a factor of 2, which is 
a convention followed in most finite element formulations. The components of the Cauchy 
stress tensor are written as; 


•[cr}- -[ Cxxj ^yyi ^zzi ^xyt ^yzt ^zxj 


(3.3) 


The matrix form of the fourth order tensor given by: 


^ J J *4 + t{^}T [CE] j 


(3.4) 


where, 


Ha} = 


Ha'} 


2 ‘a 


eq 


•A = Xn‘(ep"-‘ 


(3.5) 

(3.6) 


The array H'^H^ represents the deviatoric part of the Cauchy stress at time t and 
stands for the equivalent plastic strain at time t. For an isotropic material, the matrix 
[C^] for the 3-D case is given by: 


[C^] 



1- u 

V 

V 

0 

0 

0 


V 

1- V 

V 

0 

0 

0 

E 

V 

V 

1 - V 

0 

0 

0 

i/){l — 2i>) 

0 

0 

0 

1-21/ 

2 

0 

0 


0 

0 

0 

0 

l-2u 

2 

0 


0 

0 

0 

0 

0 

l-2u 

2 


(3.7) 


Equation (2.41) can be written in the following form owing to the symmetries of 
tAsij, tArjij and Vij: 

f 5(t{Au}^) V d^V+ [ 6{t{Asy) t{Ae} d^V 

Jty Jiy 

+ ( J(,{A,n*[r],{A,}rfV+( i(,{Aen ‘M rfV = (3.8) 

Jty J’^y 


The matrix [T] is given by: 


‘m 


*[E] 0 0 

0 ‘[E] 0 

0 0 *[E] 


(3.9) 



where, 


^XX 

t(j 

^xy 

^XZ 

<^xy 

^(7 

^yy 


^xz 

Uyz 

^CTzz 


3.2 Finite Element Formulation 


( 3 . 10 ) 


The domain is discretized into a number of elements and the incremental displacement 
field is approximated over each element by: 


(3.11) 


(Au 

t{Au} = \ ^Av I = *[$] 
tAw 

Here, the element incremental displacement vector is given by: 

t{Au}®^ = {tAu\ tAu\ tAu;\ jAu^, . . .} (3.12) 

where the quantities jAu®, jAu*, jAu;* stand for the unknown incremental displacements 
of node i in x, y and 2 : directions respectively. The matrix *[$] is defined by: 


‘[$] = 


*{^iF 


( 3 . 13 ) 


where, 


‘{$iF = fiVi,0,0,*Ar2,0,0,...} 

*{#2F = {0,*Ari,0,0,‘iV2,0,...} 

*{<&3F = {0,0,‘iVx,0,0,*Ar2,...} (3.14) 


The ^Ni, which are functions of {x,y,z), are called as shape functions. In the present 
work, 8-noded brick element with trilinear shape functions is used. The functions Wj do 
not depend explicitly on t. However, the left superscript t is used to emphasize the fact 
that the shape of element (nodal coordinates) changes with time. 

The acceleration can similarly be expressed as: 


t+At 




t+Ai^ \= 




( 3 . 15 ) 



The strain field is expressed in terms of the nodal displacements by diferentiating 
(3.11) and using the expressions (2.37) and (2.38) as: 

t{Ae} = ^[BL]t{Auy (3.16) 

t{A?7} = (3.17) 

where, 

and 

W = ■{# J,., '{# 2 },. . . .} (3.19) 

Using the equations (3.11), (3.15)^ (3.16) and (3.17), the contribution to the inte- 
gral (3.8) over a typical element with volume V® is: 

(J(t{AM}®''') V*[^] 

<J(t{A«}«^) '[Bt] of t{Au}« -t- 

5(t{A«}^^) ^{Bn] cfV^) t{Au}^+ 

<5(^{Au}«^) ^[B^f ‘{or} d = 5{t{AuY'") *+^%Fr (3.20) 

The contribution to the term is expressed in terms of the elemental external 

force vector using a standard procedure [48]. 

Since the variation in the displaceipent vector is arbitrary, the above equation can be 
written as: 

^[KYt{AuY+ W= (3.21) 

where the element mass matrix *[M]® is defined as: 

t[Mf = *[$f (3-22) 

and the elemental stiffness matrix *[KY is expressed as: 

^[KY= *[KlY+ ^[Kr^iY 


(3.18) 


(3.23) 



where 


'[KlY = (3.24) 

'[KnlY = l^/[BNf*[T]*[B^]d^V^ (3.25) 

The elemental internal force vector is; 

W= /ve (3.26) 

The element mass (‘[M]®) and stiffness Y[KY) matrices along with the element force 
vectors YfY are assembled to obtain the global equation: 

\M] + ^[K] t{AU} + Yf} = (3.27) 

where, 

‘[M] = global mass matrix 

^[K] = global stiffness matrix 

*{/} = global internal force vector at time t 

= global external force vector at i + At. 

The equation (3.27) represents a system of ordinary coupled second order linear dif- 
ferential equations. The next section describes the technique to convert it into a system 
of algebraic equations. 

3.3 Finite Difference Scheme 

The number of efficient methods for the solution of a large system of coupled ordinary 
differential equations such as those arising from a finite element formulation is limited. A 
popular technique is the direct integration of equation (3.27), where the dynamic equilib- 
rium is sought only at discrete time intervals instead of at all times t. The assumptions 
on the variations of displacement, velocity and acceleration within each time interval 
determine the accuracy, stability and cost of the solution procedure. 

Two broad classifications exist in the direct integration methods: Explicit and Implicit 
integration methods. Explicit integration methods use the equilibrium conditions at time t 
for obtaining the solution at time t+ At and hence these procedures are only conditionally 
stable. One main advantage of explicit schemes is that factorization of the effective 
stiffness matrix is not necessary. Implicit integration methods use the dynamic equilibrium 



conditions at both time t and t + At for obtaining the solution at time t + At. Since the 
implicit time integration methods are unconditionally stable, the time step used can be 
several orders of magnitudes larger, but the factorization of effective stiffness matrix is 
necessary for obtaining the solution. In this work, an implicit method called Newmark 
method is used for integration. 

The Newmark method makes use of the following assumptions [4]: 


t+At^U} 


\Lf} + [(1 - 5y{U} + 6^+^\U}] At 
'{U}+ ‘{I7}At+[(i- a)*{U}+ 


{Aty 


(3.28) . 


The parameters a and 5 can be determined so as to obtain desired accuracy and stability. 
The constant average acceleration method corresponding to a = i and | is used in 
this work, which results in an implicit unconditionally stable integration scheme. 

The Newmark constants are defined below: 


_ 1 
^(Ai)2 
1 


a{At) 



03 = At(l- (J) 

04 = 5 At (3.29) 

Substituting for from equation (3.28) and splitting as; 




(3.30) 


the expression (3.27) becomes: 

‘[ATj] ,{AC/} + (do ‘[M] ‘{U)+ ‘{/}) = ‘+“{Ci} (3.31) 

where the effective stiffness matrix and the effective force vector are given 

by: 

^[Kd] = ao *[M] + *[K] (3.32) 

^^^XFd} = *[M](oo*{f7}+ oi*{f/}+ az^M) (3.33) 

Since the dynamic equilibrium is satisfied at time t and assuming that changes in mass 
matrix in an increment is not very significant, the following equation must hold, which 



can Ue ootaxiieu noiii ine nnite elemtot discretization of the integral (2.29) written at 
time t: 

•[M] ‘{U} + ■{/} = ‘{f} (3.34) 

Using equation (3.28), this can be expressed as: 

ao\MY{U}^ ‘{/}= ‘TO (3.35) 

Decomposing the effective force vector as: 

+ t{AF4 (3.36) 

and using the expression (3.35) for *{Fd}, the following equation in t{AI7} is obtained 
from equation (3.31): 

^[Kd]t{AU}= t{AFd} (3.37) 

Using equations (3.33) and (3.35) we get: 

,{AF4= ‘+^‘TO- ‘TO= ‘{/}+ ‘[M](ai ‘TO + as ‘TO) (3-38) 

This algebraic equation system is solved to obtain the incremental displacement vec- 
tor t{Ai7}. 

Once the incremental displacement is obtained by solving equation (3.37), the velocity 
and acceleration at time t -I- At are obtained using the following expressions: 

^ aoC+At^l/} - '{U}) - ai*{U} - a2'{U} (3.39) 

^ ^{U}+ a3*{U}+ 04 (3-40) 

3.4 Determination of Stress 

The evaluation of stress components (at the Gauss points) is done by the following pro- 
cedure: 

1. The linear and non-linear strains are obtained by: 

t{Ae} = ‘TO]4A«}« (3.41) 

t{Ari} = ‘TO]4Aa}« (3.42) 

2. The Green-Lagrange strain is obtained by combining the linear and non-linear parts 
using equations (2.5), (3.1) and (3.2) in equation (2.4). 


3. FYom the Green-Lagrange strain; the Jaumann stress increment is found using the 
procedure described in section 3.4.1. 

4. Finally, the Cauchy stress increment is calculated from the Jaumann stress incre- 
ment using equations (2.2) and (2.3). This is added to the Cauchy stress at time t 
to find the Cauchy stress at time t + At. 

3.4.1 Integration of the Constitutive Equation 

Various techniques exist for the integration of the constitutive equation. A simple and 
robust technique is the Euler forward integration scheme discussed below. 

Suppose that the state of the Gauss point at time t {elastic or plastic) is known. 

• If the state at time t is elastic, 

1. Calculate the Jaumann stress increment assuming elastic behaviour. 

2. Determine the Cauchy stress increments using equations (2.2) k (2.3) and add 
it to Cauchy stress at time t, *{cr}, to get the the Cauchy stress at time t-\-At, 
viz., 

3. Detei'mine and using equation (2.10). 

4. If < ‘cTy, then the material point is still in elastic state. Return. 

Else, a transition from elastic to plastic has ocurred. Calculate: 

t+At _ 

(3.43) 

and change the state to plastic. The sub-incrementation method is followed 
and the Jaumann stress increment is updated after each sub-increment by the 
change in stress components corresponding to the elasto-plastic strain sub- 
increment. Thus, the Jaumann stress increment is given by the relation: 

t{a At} = [C^]{1 - Ratio) t{Ae} + {Ratio/n) t{Ae} (3.44) 

Z=1 

where the matrix corresponds to (z - l)^'* updated state. 

• If the state at time t is plastic, check for unloading as described in section (3.4.2). 
If there is no unloading, the sub-incrementation method described in (4) above is 
applied with Ratio set to unity. 


3.4.2 Unloading Scheme 


The phenomenon of local unloading is to be incorporated to reproduce more closely the 
elasto-plastic response of a structure. The unloading criterion defined by Chakrabarty [8]; 


‘{nrdAe}<0 


(3.45) 


where, represents the unit outward normal to the yield surface at the current stress 
point V in a 9-D stress space, is implemented in this work. Neutral loading and postive 
loading are represented by replacing the “<”symbol in (3.45) by “=” and “>” symbols 
respectively. 

When unloading is detected, the yield stress at that Gauss point is changed to the 
equivalent stress at the Gauss point at time t {*ay = the state tag of the Gauss 

point is changed from plastic to elastic and the elastic constitutive equation is applied to 
calculate the incremental stress. 


3.5 Skyline Scheme of Global Assembly 

The three-dimensional dynamic elasto-plastic contact problem is computationally inten- 
sive one and therefore every effort should be made to reduce the storage and the compu- 
tational time required for solving the problem. Skyline assembly scheme [4] is one such 
technique, which significantly reduces, the storage required for the global stiffness ma- 
trix ^[K]. In this storage scheme, only the elements below the skyline of the matrix ^[K] 
are stored in a one-dimensional array *{A}. However, along with the storage scheme a 
specific procedure for addressing the elements of *[K] in ^{A] is needed. Thus, before 
proceeding with the assembly of the element stiffness matrices it is necessary to establish 
the address of the global stiffness matrix elements in one-dimensional array 
The element pattern of a typical global stiffness matrix is shown below: 

kn ki2 0 ki4 0 

ki2 ^22 ^23 0 0 

^[K] = 0 A;23 kzs k^i 0 

kii 0 k$4 kii /s45 

0 0 0 A;45 ks^ 

Since the matrix ‘[isT] is symmetric, only the elements above the diagonal and the diagonal 
elements need to be stored. Defining by m the row number of the first non-zero element 



in column i, the variables rrii = represent the skyline of the matrix, and the 

value {rrii+i — mi) represent the column height. (The zero elements below the skyline are 
stored, since they may become non-zero during the solution process). The column heights 
are determined from the connectivity arrays of the elements. Once the column heights 
of the global stiffness matrix are known, we can store all the elements below the skyline 
of *'{K\ as a one-dimensional array in *{A}; i.e., the active elements [elevaents below the 
skyline) of ^[K] are stored consecutively in 

{^Hj ^22) ^12) kzZi ^23) ^44) ^34) Oj ^14; ^55) 1^45^ 

In addition to ‘{A}, we also define an array {MAXA}, which stores the addresses of 
diagonal elements of ‘[iiT] in *{A}; i.e., the adress of the diagonal element of *[iif], ku, 
in array *{A} is given by MAXAi. Thus, 

{MAXA}^=:{1,2,4,6,10, 12} 

Note that, M AXAi is equal to the sum of column heights upto the (i- 1)*'^ column plus 1. 
Hence the number of active elements in the column of *[K] is equal to MAXAi+i — 
MAXAi, and the element addresses are MAXA^, {MAXAi + 1), • • •, (MAXAj+i - 1). 
It follows that using the above storage scheme of ‘[iiT] in *{A}, together with the address 
array {MAX A) as defined above, each element of ^[K\ in *{A} can be addressed easily. 

3.6 Static Condensation Scheme 

The contact analysis is done in an iterative manner to determine the nodes which are 
out of contact, in sticking contact or in slipping contact. Solving the complete system of 
equations in each iteration is highly inefficient in terms of time and computer resources, 
especially when the contact region is small compared to the whole system. In order to 
reduce the solution time, it is desirable to condense the effective stiffness matrix and the 
effective force vector to the size of the degrees of freedom to which the contact boundary 
conditions are to be applied. 

The condensation procedure starts by dividing the nodes into two categories: the nodes 
at which the contact boundary conditions are to be applied (referred as the nodes of type 1 ) 
and the remaining nodes (referred as the nodes of type 5). To facilitate the condensation, 
the nodes are renumbered so that the nodes of type 1 are numbered first (called internal 



node numbers, while the origmal node -numbers are called external node numbers). Then 
the essential boundary conditions (i.e., the displacement boundary conditions) are applied 
to equation (3.37). 

Then these equations are partitioned as follows;^ 


’[Jrii) [Kn] 

< 

1 _ , 

' {AF.} ( 

^K 2 l] [K 22 ] _ 


i {A!/,} J 

[ {AF,} 


(3.46) 


where {At/j} denotes the incremental displacement vector corresponding to the nodes 
of type 1 and {AL/^} denotes the incremental displacement vector corresponding to the . 
nodes of type 2. Equation (3.46) can be seperated as follows: 


[■^u]{AC/i}+ [ii:i2]{A172} = {AFi} (3.47) 

[K 21 ] {AC/i} + [K22]{MJ2} = {AE 2 } . (3.48) 


Solving for from equation (3.48), we get: 


{At/ 2 } = [K 22 V ({AF 2 } - [K 21 ] {AU,}) (3.49) 


Substitution of this expression for {AC/ 2 } in equation (3.47) and rearrangement of the 
resulting equation leads to the following condensed set of dynamic equilibrium equations: 


[K]{AU} = {AF} (3.50) 

where, 

[K] = [Kn]-[K,2][A] (3.51) 

{AC/} = {AUi} (3.52) 

{AF} = {AFJ - [F 12 ] {5} (3.53) 

[A] = [F22]-'[^2i] (3.54) 

{B} = [F22]~HAF2} (3.55) 


The nuinber of equations in the condensed set (3.50) is equal to the number of degrees 
of freedom of type 1. The coefficient matrix [K] and the right side vector {AF 2 } are 
evaluated from equations (3.51 - 3.55) using the partitioned matrices and vectors of 
equation (3.46). Note that while evaluating [A] and {B}, it is not necessary to invert the 


Tn this section, the superscript /subscript denoting time have been omitted for the sake of convenience 



matrix [^^ 22 ] • Instead one can solve the following equations by the Gauss elimination 

method: 

[i^22] [-A] = [K 21 ] (3.56) 

[K-izliB] = {AF 2 } (3.57) 

To reduce the time further, one can store the upper triangular form of [K 22 ]. In 
this way, one can solve for the columns of [A] by performing Gauss elimination and 
back substitution operations on the corresponding columns of [K 21 ]. The vector {B} is 
obtained in similar fashion by performing the Gauss elimination and back substitution on 
the vector 

3.7 Modified Newton-Raphson Scheme 

The solution to equation (3.37) represents only an approximate solution to the governing 
equation (2.40), because of the linearization and simplifications employed in deriving the 
equation (3.37). The error obtained should to be minmized and there are various iterative 
techniques for this. One of the simplest and very effective technique is the Modified 
Newton- Rap hfi on algorithm, which offers fast convergence with less computation [4]. 

This tecbniciue can elegantly be represented as follows: 

Solve: 

^ [K,] (3.58) 

where, 

(3.59) 

= t{AF4 (3.60) 

till the convergence criterion: 

l|.{AFu 

is satisfied. 'I'he vector is called the unbalance force vector. 



3.8 Divergence Handling Procedures 

The modified Newton Raphson method diverges in some cases; i.e., the ratio of unbalance 
force to the external force, instead of reducing, increases in each iteration. Also, in 



some other cases, the rate of convergence is not fast enough, requiring a large number 
of iterations for reducing the unbalance force to acceptable limits. Two of the methods, 
which are simple and fairly effective for handling divergence or accelerating the rate of 
convergence are Under-relaxation method and Line search method [7] . These methods 
are generally applied by modifying the displacement vector However, in the 

contact analysis scenario, it is not acceptable, since the contact status may get changed 
due to modification in the displacement vector. Therefore, in this work, these methods are 
applied by modifying the unbalance force vector and repeating the solution 

process when divergence is detected. 

In this work, the Newton Raphson iterations are said to be diverging when; 

(3.62) 




is satisfied. The value of the factor /3 is taken as 0.9. 

The under-relaxation and line search methods are discussed below: 

Under-relaxation method: If divergence is detected in the i*^ iteration, the unbalance 
force v(K;t<)r ({A/'u}h"’) is scaled by a factor a repeatedly, till the above criterion is 
not satisfied: 


= a t{AT;}(*-') (3.63) 

If the under- relaxation is not able to find a value of unbalance force for which 
the above criterion is not satisfied, the solution run is terminated with a warning 
messag('. 


Line search method: In this method, the basic aim is to search for the point, which 
corresponds to minimum unbalance, along a direction specified by the unbalance 
force vector < { However, it is not very efficient to perform a full line search, 

since, in that case the number of solution runs required may become prohibitively 
high. Th(i approach used is to repeat the solutions for a fixed number of values of 
factor otj (defined below), within a preselected range. 

(3.64) 

where, ^min < otj < ^max 


The value of aj corresponding to the minimum unbalance force is chosen as the 
best point and the solution corresponding to this (Xj is accepted as the result of line 


search method. 



Since in both the above methods, the complete iteration is repeated, including a 
number of contact iterations, it is very costly to do either under-relaxation or line search, 
in terms of both computational time and resources. Therefore, it is preferable to use 
these methods only when some difficulties in convergence are observed. Under- relaxation 
is cheaper compared to line search, while line search is more efficient. 

3.9 Numerical Aspects 

The selection of parameters controlling the numerical scheme in dynamic non-linear anal- 
ysis is of imporatance not only for efficiently obtaining the solution, but also for the 
accuracy of the solution. The following sections discuss various aspects related to the 
numerical scheme used. 

3.9.1 Choice of Time-step 

A procedure for the selection of time step is given in Bathe [4], which may be used to select 
an initial time step and mesh. A time step and mesh convergence study should then be 
performed to arrive at the proper values. It is important to note that an unconditionally 
stable method does not provide a solution close to the exact solution for an arbitrary time 
step. Phenomena such as amplitude reduction and period elongation may occur due to 
an improper choice of time step. Hence the proper way to perform a dynamic non-linear 
analysis is to study the effect of time step and mesh size on the solution and to search 
appropriate values. 

The choice of a time step is based on the type of problem considered. In general, 
dynamic problems can be categorized as: 1. Strucural dynamics problems and 2. Wave 
propagation problems. 

Structural Dynamics Problem: The step by step procedure for modeling a strucutu- 
ral dynamics problem is given below. 

1. Identify the frequencies involved in loading, using a Fourier analysis if nec- 
essary. These frequencies may change with time. Let the highest significant 
frequency contained in the loading be Wu- 



2. Choose a finite element mesh that can accurately represent the static response 
and accurately represent all frequencies upto about Wco = where Wa, is 
called as the optimum critical frequency. 

3. Perform the direct integration analysis. The time step At for this solution 
should be equal to about where T^, = 

Wave Propagation Problems: Assuming that the critical wavelength to be calculated 
coirectly is then the time taken by the wave to travel past a point is given by: 

^ (3.65) 

where c is the wave speed. Assuming that n time steps are necessary to represent 
the travel of the wave, the time step should be; 

At = ^ (3.66) 

and the effective length? of a finite element should be: 

Le = cAt (3.67) 

The above choice of time step and corresponding effective length represents the 
complete wave travel accurately. But, care has to be taken to ensure that At < 
where 'i’„ is the smallest time period of the problem. 

It is found that the most accurate solution is obtained by integrating with a time 
step equal to the above limits (denoted by Ata) and the solution is less accurate, when 
a smaller time step is employed. This deterioration in the accuracy of the predicted 
solution when At is smaller than Ata- is more pronounced, when a relatively coarse 
spatial discretization is used. 

3.9.2 Deformation Dependent Loading 

Deformation dependent loading is taken care of by evaluating the consistent force vec- 
tor at time t -h At for the calculation of the unbalance force vector (3.59). An 

alternate procedure given by Bathe et al [11] incorporates deformation dependency into 
the formulation of the governing equation, but this results in an asymmetric stiffness ma- 
trix, the solution of which is computationally expensive. The rnethod used in this work 


^Effective length is the smallest distance between any two nodes in the mesh used 



may take a slightly more number of iterations for convergence, but this loss is oifset by 
the gains acquired by retaining a symmetic stiffness matrix. 

3.9.3 Stress Updation 

There exist two methods for updating the stress at a Gauss point: 1. Iterative updation 
and 2. Incremental updation (Crisfield [7]). 

Iterative Updation: The iterative displacement is used to compute the iterative strain, 
which is then used in the constitutive equation to update the stress at iteration (z-1) 
to stress at iteration (i). 

Incremental Updation: All the iterative displacements upto iteration (z) are added to 
find fh(^ cumulative incremental displacement of the current increment. This is used 
to C()mj)ut,e the incremental strain, which is used in the constitutive equation to 
updat e the stress from last equilibrium point (end of the previous increment) to the 
state after iteration (i). 

The first method may lead to spurious unloading, since the unbalance force vector 
calculated in some iterations may be in different directions compared to that of the initial 
force vector. This can be avoided by using the incremental updation procedure. 



Chapter 4 


Contact Formulation 


The dyiianiic-, laige defoi rnatioii, elasto-plastic, updated Lagrangian finite element formu- 
lation developed in the last two chapters, can be used for contact analysis with appropriate 
modifications. 1 he most important modification is the development of an additional set of 
finite element ecjuations (involving the contact stiffness matrix) relating the unknown con- 
tact forces and displacements. These equations, which consist of the kinematic constraints 
and the con(,act force expressions based on the discretization of the contact boundary, are 
developed using the Lagrange multiplier method. These developments are described in the 
following sections. Finally the algorithm for the analysis of dynamic, large deformation, 
elasto-plastic contact including the effects of friction, is described. 

4.1 Contact Constraints 

A typical two body contact system is shown in figure (4.1). The bodies occupy the domains 
i+A(^5 and at time t -+■ At. The boundaries of and are denoted by 

t+Atsi volumes by and respectively. At any 

time t + At, the boundary of contact body n (n = 1,2) can be partitioned as; 

t+At^n ^ t+Atgn jj t+Atgn y t+Atgn 


where, 

— surffice on which displacements are specified, 

= surface on which tractions are specified, 

‘+^‘S'c = surface on which contact occurs. Assuming the boundary to be 




Figure 4.1: Two body contact system 

smooth, outward unit normal vector at point^ f+A^n boundary is denoted 

by Other two orthonormal boundary vectors (in tangential directions) are 

and (see Figure (4.2)). 

Suppose that two boundary points and *+'^^2 ^^j-g contact with each other at 

time t + At and the contact traction at is denoted as Then by Newton’s 


third law of motion, we have: 


t+Atfj^l t+At'jp2 


where can be expressed as: 

3 

t+Atfjrn _ t+At^n t+Atj^ (4 3 ) 

* 1=1 

where, denotes the component of along 

If the two contact bodies are not welded together, then the normal component of 
traction can not be tensile. Thus we have: 


t+At^n < 0 

Further, from the Coulomb friction law we have: 


/(t+Atin)2+ (t+Atfn)2 < 


^Vectors are denoted by a letter with underbar 



Hitting Body (Body 1) 



Figure 4.2: Points in contact and associated normal and tangent vectors 

In equation (4.5), the “<” symbol represents the sticking friction and “=” represents the 
slipping condition. Further, the resultant tangential (friction) force is in the direction 
opposite to the relative velocity of slipping. The equations (4.4) and (4.5) represent 
traction constraints. 

Physical considerations require that one body can not penetrate the other. Thus, for 
the points in contact, pentration must be zero. This condition can be mathematically 
stated as: 

t+At^2 ^ Q (4.6) 

where represents the penetration between the two points in contact. In addition 

to (4.6), we need to have two more conditions to enforce the sticking friction, which 
ensure that there are no relative displacements (slip) between the contacting points in the 
tangential directions. These conditions can be stated as: 



= (t+A^2 _ t+At^l). t+At^2 ^ Q 

(4.7) 


_ («+A^2 _ t+At^l)_ t+At^ ^ Q 

(4.8) 


The conditions (4.6), (4.7) and (4.8) are called as the kinematic constraints. 



4.2 Contact Force Expressions 


As the contacting bodies aie discretized for the purpose of developing finite element 
equations, the contact boundary gets automatically devided into surface elements, called 
as contact segments. The contact force formulation in this work employs a node-to-surface 
interface model, since this work deals with a large deformation problem including slipping 
at the contact interface. Normally in a contact formulation, contact is considered only at 
discrete nodes. Here, it is assumed that a node of the hitting body (body 1) comes into 
contact with a segment of the target body (body 2). The first step in the development of 
the contact force expression is the determination of the unit normal and tangent vectors 
at a point on the target segment. This is described in the next paragraph. 

Geometry of a target segment can be described using 2-D shape functions. Thus, 

= (4.9) 

3=1 

where, are the shape functions, denotes the nodal position vector and N 

is the number of nodes per target segment. The side of a segment on which contact can 
occur is referred to as the positive side of the segment and the other side as the negative 
side of the segment. The local numbering of the nodes on a contact segment should be 
counter clockwise when we face the positive side of the segment. This is to get the correct 
sense for the outward normal vector at given point. A normal vector at point rj) 

on the positive side of the segment is defined as: 


t-hAtjsj-2 _ t-hAtjY^ ^ t-i-AtjY^ 

(4.10) 

where, and are the two tangent vectors at given by: 

t+Atjy2 ^ 

(4.11) 


(4.12) 

The unit normal and tangent vectors are defined as: 


t+At^2 _ t+AtjY2^|t+Atjy-2| 

(4.13) 

HAtj^2 ^ t+At^y|t+Ai^| 

(4.14) 


(4.15) 

Since the contact is considered only at discrete hitting nodes, we can deal only with 

point forces at hitting nodes rather than with contact tractions. 

At time t + At, let 



f+At{^ 2 } tje an array of Cartesian components of point force vector at a point of 

the target segment at which the hitting node makes contact. By Newton’s third law, the 
force vector on the hitting node will be — 

Let } and {(5^u } be the virtual displacement vectors at time t + At at 

these points (on the target and hitting bodies respectively), when they come in contact. 
Then the virtual work at hitting node can be expressed as:^ 

6 wc= (5 (4.16) 

Note that the elements of namely are the virtual 

nodal displacements of the hitting node. The elements of can be expressed in 

terms of the virtual nodal displacements of the target segment 
i = 1, N) using the 2-D shape unctions ($j, i = 1, N). Thus; 

^ (4.17) 

where, 

-10 0 0 0 #2 0 0 . . . 0 0 

[Qc] = 0 -1 0 0 $1 0 0 #2 0 ... 0 0 (4.18) 

0 0 -1 0 0 0 0 $2 . • • 0 0 $JV 

and 

(4.19) 

4.2.1 Sticking Friction Condition 

For the sticking friction, the vector can be expressed as: 

t+Atj2 
t+Atj2 

or, 

t+At^p2y _ t+At^j^2y+At^^2'^ ( 4 . 21 ) 

where, represents the Cartesian component of the unit vector *'*' and /j 

are the components, with respect to of the contact force vector *'*' *£ . Combining 

^In this section, the right superscript denoting the iteration number has been omitted for convenience 



i+Atr Tr2\ _ 


t+Atfq^ t+At^^ 
t+Atjv2^ t+At^2^ ‘+^*A/|2 



the equations (4 .16), (-4-17) and (4.21), we get the following expression for the virtual work 
at the hitting node: 


5wc= *+^‘{(5 

(4.22) 

where the vector is given by; 



(4.23) 


This is the contact force expression for a sticking node (hitting node, which is under 
sticking friction condition). 

4.2.2 Slipping Condition 

The contact forc.e expression for a slipping node (hitting node, which is slipping on the 
target segment ) is developed below. 

The vector ^ hi equation (4.23) can be split as a contribution from the normal 

contact, forci's and tangential contact forces. Thus, 

( 4 - 24 ) 


where. 


HAt 


C»} =.10of ‘+“/f 


(4.25) 


and 

< .. = [Qf (‘+“|jv|} ‘+"/| + ‘+“{JV|} ‘+“/|) (4.26) 

Here the arrays ‘ = ^>^>3 contain the Cartesian components of the vec- 

tor 

The three components of contact forces are not independent for a slipping node. Prom 
Coulomb’s friction law, we have the relation; 

^(t+At/2)2 + (t+Atfiy = (4.27) 

where p, is the coefficient of friction, for a slipping node. Therefore, the tangential com- 
ponents of the contact forces can be written as; 

t+Atf2 _ _ *+Aty2 ^ (.Qg 0 

t+Aty2 _ _ t+Atj:2 ^ gjjj 0 


(4.28) 

(4.29) 



where, 9 is the angle between the resultant tangential (friction) force and the vector 
and since the noinial contact force is always compressive, ” sign is used in the right 
hand side of both the above equations. 

Calculation of the angle 9: 

The friction force in the case of slipping is in a direction opposite to the relative 
velocity of slii)j)ing, which is same as the difference in incremental displacement of the 
target point and the hitting node. Therefore, the direction of the friction force on the 
target body is given by: — = — ([QJ t{Auc})j, where the right subscript t 

indicates the tangential component. Define: 

= -‘+^‘{iV2T[Qc]t{Aac} (4.30) 

,A53 = - [QJ ,{A«4 (4.31) 


Thus, the angle 0 is given by: 


9 = tan ^ 


Auz \ 

tAu2y 


(4.32) 


However, to avoid the difficulty in handling infinity, the computer implementation uses 
cosine function instead of tangent function for calculating 9. 


Using ecpiations (4.28) and (4.29) in equation (4.26) and substituting equations (4.25) 
and (4.26) in ecjnation (4.24) we get: 

} = [ClY cos 9- la sin 9) (4-33) 


4.3 Kinematic Constraints in Nodal Form 

The kinematic constraints should be written separately for the sticking nodes and for 
slipping nodes. 

4.3.1 Sticking Friction Condition 

To express the kinematic constraints in terms of the nodal displacements, the equations 
(4.6), (4.7) and (4.8) are combined in matrix form for the i^^ Newton-Raphson iteration 
as: 

t+At{p}(i3 ^ t+At|2;2 _ 


(4.34) 



where, 




‘{p}“> = (4.35) 

We have, 

= '+^>{12 _ ,jl}|i-l) ^ jgj- 

where, the seconci term on the right hand side is an array of Cartesian components of 
the difference of the displacement vectors of the target and hitting bodies at the contact 
node, for the current iteration. For the first iteration: 




(4.37) 


Analogns to equation (4.17), this can be written as: 

ifAn'-" - = [Q,] t{Aue}(< 


whoni, the vector ({Aiq.}^’^ given by: 


:{An<;} 


a? 


{,Ak"‘>, .AC™, Auf". ,A„J'''. ,hwf\ 


,1(0 


ih) a„,2(0 


„2(d 


„2(*) 


(4.38) 


lAffiJ,'*’} (4.39) 


contains tlui nodal values of displacement at the hitting node and the nodes of the target 
segment for the current iteration. Recognizing the first term of the right hand side of the 
equation (4.3G) as corresponding to (penetration for the previous iteration), 

equation (4.30) can be written as: 




-{p}h)= 

For the first iteration, 

«+At{p}(i-l) _ HAt|pj.(0) _ 

Tliereforo, the kinematic constraints for sticking node is: 

+ ‘+“[iV^]'’[Qd .{A«J<‘> = {0} 


j) 


i+/ 


(4.40) 


(4.41) 


(4.42) 


The first term of the right hand side of equation (4.36), on similar development gives: 


t+A£|-pj(i-l) _ t+A£ 




(4.43) 


where, is the array of coordinates of the hitting node and of the nodes of the 

target segment updated upto the previous iteration. 



4.3.2 Slipping Condition 

The development of the nodal constraints for the slipping nodes is similar to the sticking 


friction condition, except that we have only one kinematic constraint, viz., the equa- 
tion (4.6). 

We can write the matrix form of the kinematic constraints for the 

iteration as: 

Newton-Raphson 

t+A^P ^ t+At^jy2jrt+Atl^2 _ 

Using the equations (4.35 - 4.39) in (4.44) we get: 

(4.44) 

For the first iteration. 

(4.45) 

t+A^(i-l) ^ t+A^(0) ^ 

Therefore, the kinematic constraint for a slipping node is: 

(4.46) 

Similar to the sticking condition, we have: 

(4.47) 


(4.48) 


where, contains the coordinates of the hitting node and the nodes of target 

segment update<I upto the previous iteration. 

4.4 Lagrange Multiplier Method 

In Lagrange Multiplier Method, the contact forces are considered as primary unknowns 
and the kinematic constraints are enforced exactly. We have the contact force expression 
for a sticking node given by equation (4.23) and for slipping node by equation (4.33). The 
kinematic constraint for sticking node is given by equation (4.42): 

[JV=] (OJ «{ (4-49) 

and for a slipping node is given by equation (4.47): 

<+A.{^2jr jQj .{AjiJW = 


(4.50) 



(4,51) 


Combining the equations (4.23) and (4.49) to get for a sticking node, 


"“fcj ' 


^ = < 


0 

1 J 




where, 

Combining the equations (4.33) and (4.50) we get for a slipping node, 


0 - ‘ 


[ ■+A</2(‘) 1 1 


1 

o 

< 


[ J ” 1 



where, 


(4.52) 


(4.53) 


"“{'ft} = Wcf ("“{Alf}- "“{Af|}itcos9- ‘+A‘{iV|},usm«) (4.55) 


Assembling these equations over all the potential contact nodes, we get the following 
global equation: 

0 - 

_ - 0 

Note that, while the whole global coefficient matrix is known from the geometry, only 
a part of the right side vector, namely is known from the geometry. The 

other part, is unknown. This vector is eliminated by combining this equation 

with equation (3.58) of chapter 3. Before this, the set of equation (3.58) is condensed 
to retain only those degrees of freedoms, which are associated with the contact nodes. 
This condensed and combined set is solved to find the nodal displacements of the contact 
nodes and the contact forces. Then the equation (3.49) is used to find the displacements 
of nodes of type 2. 



r 4+At|-pj.(j-l) 

. .{At/J(') J 

1 


(4.56) 


4.5 Algorithm used for Dynamic Large Deformation 
Elasto-plastic Contact Problem 

For solving a typical dynamic, large deformation, elasto-plastic, contact problem, the 
following steps are used: 



1. Finding the potential contact nodes for which the hitting and target nodes are clos 
than a prescribed length. 

2. Renumbering the nodes: Nodes of the hitting and target nodes are renumbered su 
that the contact nodes are numbered first. This is carried out to facilitate the stat 
condensation of the stiffness matrix which helps in reducing the computational tim 

3. Coefficient matrix: Based on the current geometry and state of stress, condenst 
form of effective stiffiness matrix and force vector is formed. 

4. Begin the Newton-Raphson iteration. 

5. Contact search: Search for the target segment corresponding to each hitting nod 
Master-slave algorithm is used here. 

6. Contact iterations: 

• Initially all the potential nodes are assumed to be in sticking friction. Fori 
the contact stiffness matrix and right side vector. Combine with the condense 
form of the effective stiffness matrix and the effective force vector. Solve th 
system of equations. 

• Find the “out of contact” nodes, for which the normal component of the contac 
reaction is tensile (refer equation (4.4)). The contact reactions at these node 
are forced to be zero during the subsequent contact iterations. 

• After removing all the “out of contact” nodes, the nodes which are slipping ar 
determined. We know that a node is slipping, if the contact reactions at tha 
node violate equation (4.5), expressed in a nodal form. 

• Repeat the contact iterations till the correct direction of friction force is oh 
tained for all the slipping nodes and no node changes the contact status. 

7. Find the displacements of the non-contact nodes (type 2 nodes) from the displace 
raents of type 1 nodes. 

8. Updating: Update the stresses, strains and contact forces. 

9. Convergence: Find out the unbalance force and check for convergence. If th 
Newton-Raphson iterations have converged, store the results at the end of the incr€ 


ment and start the next time increment. If the convergence criterion is not satisfied, 
repeat the Newton-Raphson iteration and continue till convergence. 



Chapter 5 


Results and Discussion 


A finite element code for the analysis of 3-D, dynamic, elasto-plastic contact problem 
including the effect of friction is developed based on the formulations presented in the 
chapters 2, 3 and 4. The stiffness matrix, mass matrix and force vectors are evaluated 
using numerical integration, while Newmark’s scheme is used to update the velocity and 
acceleration after each time-step. The code is developed for 8-noded brick element and 
2-point Gauss-Legendre numerical integration is used in each direction. This chapter 
discusses the validation of the code and a study of an additional problem to show the 
application of the code to impact problems. 

5.1 Validation 

The validation of the code is necessary to prove the applicability of the program to the 
class of the problems of interest. Here, the validation is done by comparing the results 
obtained with the results reported in the published literature. 

The code developed for dynamic large deformation elasto-plastic contact analysis is 
validated by using the following problems; 

1. Elasto-plastic dynamic response of a simply supported square plate. 

2. Oblique impact of two infinite blocks. 


3. Impact of two cantilever beams. 



5.1-.1 Elasto-plastic Dynamic Response of a Simply Supported 
Square Plate 

This is a single body analysis of an elasto-plastic, sinaply supported square plate subjected 
to a lateral step load of 300 psi, reported by Gendy and Saleeb [49], The material is 
assumed to be elastic-perfectly plastic. The geometric and material properties of the plate 
as given in the reference [49] are given below: 

• Length of side = 10 inches 

• Thickness = 0.5 inches 

• Young’s modulus = 1 x 10^ psi 

• Poisson’s ratio = 0.3 

• Yield stress = 3 x 10^ psi 

• Density = 2.588 x lO'^ Ib-secVin^ 

Due to the symmetry, only one quater of the plate was modelled. The finite element 
iriodel consist of 10 elements along the side (half length of the plate) and 4 elements along 
thie thickness. The timestep used is 22.3 x 10”® seconds, which is approximately equal to 
( 1 /A&y^ of the fundamental period of the linear elastic plate. 

Since, the material model used in this work is based on a power law type hardening, 
t can not be used for an elastic-perfectly plastic material. However, the response for 
h.e elastic-perfectly plastic material can be obtained as the limiting case when iT 0. 
'ence for different cases are considered, where K is reduced progressively: elastic material, 
asto-plastic material with iT = 1 x 10® psi, n = 1, elasto-plastic material with K = 

< 10® psi, n = 1 and elasto-plastic material with iT = 1 x 10^ psi, n = 1. For K smaller 
m 1 X 10^ psi, not much change is observed in the response. Hence, this last case can 
considered as the limiting case. The displacement at the center of the plate given in 
-ence [49] and those obtained using the present code are shown in figures (5.1) and 
respectively. An excellent agreement in trends of the predicted time responses is 
- J. However, the displacements obtained in this work are slightly lesser than those 
u f-TG]. This is mainly due to two reasons: 



1. numerical inaccuracies associated with approximating an elastic-perfectly plastic 
material using a hardening type material model and 

2. different types of finite elements used (Gendy and Saleeb [49] used higher order shell 
elements, while linear 8-noded brick element is used in this work). 

Figures (5.3) and (5.4) give the variation of maximum equivalent stress and maximum 
equivalent plastic strain with time. It is observed that the stress increases slightly even 
after reaching the yield point, which is not the case for an elastic-perfectly plastic material. 



Figure 5.1: Displacement at the center of the plate (Ref. [49]) 



Figure 5.2: Displacement at the center of the plate (Present work) 
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Figure 5.4: Maximum equivalent plastic strain in the plate (Present work) 



5.1.2 Oblique Impact of Two Infinite Blocks 



Figure 5.5: Oblique impact of two elastic blocks 

The second problem used for validation is the oblique impact of two infinite elastic 
blocks (plane strain problem) shown in figure 5.5. This problem is reported in refer- 
ence [43]. The top block is given a uniform initial velocity, Vq = (—10,0, —10) units (in 
the coordinate system shown in figure). The geometric and material properties are; 

• Hitting Block 

— length = 5 units 

— height = 4 units 

• Target Block 

— length = 16 units 

— height = 6 units 

• Young’s Modulus {E) = 1000 units, for both the blocks. 

• Poisson’s ratio {v) = 0.0, for both the blocks. 

• Density (p) =0.1 units, for both the blocks. 



Armero and Petocz- [43] have solved this problem as a plane strain problem. Bhat [50] 
has solved the same problem, only for the frictionless case, using Jaumann stress measure 
and Green-Lagrange strain measure. Since the Poisson’s ratio is zero, the plane strain and 
plane stress formulations will give the same results. In this work, the problem is solved 
with only one element in the F-direction. To overcome the numerical difficulty arising 
due to low stiffness in F-direction, the displacement along this direction is fixed at all 
nodes for both the bodies. 

The following are the discretization and timestep considered: 

• Discretization; 

— Number of elements in the hitting block = 20 (5 x 1 x 4) 

- Number of elements in the target block = 96 (16 x 1 x 6) 

- Number of d.o.f in the hitting block = 180 

- Number of d.o.f in the target block = 714 

• Time-step = 0.01 units 

There are two different cases. In the first case, there is no friction at contact interface and 
in the second case friction coefficient (yu) between the hitting and target bodies is 0.4. 

The X and Z displacements at point P are plotted against the time. The F direction 
used in reference [43] is the .2'-direction in this work. Both the displacements and time 
are measured from the instant of contact between the two blocks. The results reported 
in [43] are given in figures (5.6) and (5.8) and those obtained from the present code are 
shown in figures (5.7) & (5.9). The results are found to be in good agreement. 
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Figure 5.7: X-displacement of point P (present work) 




Figure 5.8: Y-displacement of point P (Ref. [43]) 



Figure 5.9: Z-displacement of point P (Present work) 






Figure 5.10: Impact of two cantilever beams 

5.1.3 Impact of two Cantilever Beams 

This is a large (it'forniation, elastic, plane stress problem given in reference [35]. Both the 
cantilever Ixuims arc ichmtical and fixed at different ends as shown in figure (5.10). The 
geoinelric. and material properties are given below: 

• Length of cantilever = 20 units 

• Height of cant ilever = 1 unit 

• Thickiuiss of cantilever = 0.01 units 

• Density = 4 units 

• Young’s modulus = 4000 units 

• Poisson's ratio = 0.3 

There are two cases of external loading P, and P„ each with a magnitude 0.15 units and 

three cases of friction (/r — 0.0, 0.2, 0.4) for each case of loading. 

Sung aiKi Kwak [3,5) have done the problem using plane stress formulation. To do 
this problen. using three dimensional problem, we have to prevent the twisting of the 
cantilever, produced due to the errors associated with numerical computation. This ob- 
jective is achieved by modelling only half of the cantilevers and applying the symmetry 

boundary conditions at the mid-plane. 

The discretization and time-step size axe given below: 

• Number of elements for each cantilever — 80 (40 x 1 x 2) 



• Number of dof for each cantilever = 738 


• Time-step size = 0.015 units 

The Horizontal and vertical displacements of the end tip of the upper cantilver reported 
in Sung and Kwak [35] are given in figures (5.11) and (5.13) and those from obtained in 
this work are given in figures (5.12) and (5.14). It is found that, for the frictionless 
case, the horizontal displacements axe closely matching till time = 10 units. Beyond this 
time, the results from this work are substantially higher than those reported by Sung 
and Kwak [35]. The agreement is better for the vertical displacements, which show close 
matching with the reported results till time = 14 units for case PII and till divergence 
(time = 11.25) for case PI. The difference in the results in the later time incerements may 
be due to the numerical errors accumulated over the increments and due to the differences 
created by plane stress and 3-D formulations. 

For the cases where there is friction, convergence difficulties are observed, even after 
using the divergence handling algorithms. The analysis could be completed only for the 
loadcase PI, with coefficient of friction fx = 0.2. The results obtained for this case also 
are shown in figures (5.12) and (5.14). The results show good agreement for horizontal 
displacement till time = 8 units and for vertical displacements till time = 14 units, with 
those reported by Sung and Kwak [35]. 




Figure 5.11: Horizontal displacement at the end tip of upper cantilever (Ref. [35]) 



Figure 5.12: Horizontal displacement at the end tip of upper cantilever (Present work) 
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Figure 5.13: Vertical displacement at the end tip of upper cantilever (Ref. [35]) 



Figure 5.14: Vertical displacement at the end tip of upper cantilever (Present work) 




1 




1 


1 


Figure 5.15: Impact of a sphere on plate 

5.1.4 Elasto-plastic Impact of a Sphere on a Plate with Friction 

To the demonstrate the applicability of the program developed, elasto-plastic impact of 
a sphere on a plate is analyzed (See figure (5.15)). Two different cases, viz., frictionless 
contact and contact with friction interphase (/x = 0.2) are considered. The sphere is given 
an initial velocity, V = (10,0,-10), in the coordinate system shown. 

The geometric details of the sphere and plate are shown in figure (5.15). The material 
properties assumed for both the bodies are as follows: 

• Young’s modulus (E) = 4 x 10® units 

• Poisson’s ratio {u) = 0.3 

• Density (p) = 4 units 

• Yield Strength ((cry)o) = 5000 units 

• Hardening coefficient (K) = 40000 units 


Hardening exponent (n) =0.4 



Figure 5.16: Finite element plot of sphere (Full view) 


The plate is fixed along the edges. Time step used for the analysis is 0.002 units. The 
mesh details of the plate and sphere are as follows: 

• Number of elements for sphere = 672 (Refer figures (5.16) and (5.17)) 

• Number of elements for plate = 1024 (16 x 16 x 4) 

• Number of d.o.f. for sphere = 2319 

• Number of d.o.f. for plate = 4335 

The plots of the configurations of sphere and plate at different times are given in 
figures (5.18 - 5.21) for frictionless case and in figures (5.22 - 5.25) for the frictional case. 
The X and Z displacements of the bottom node and the centre of the sphere for both the 
frictionless and the frictional cases are shown in figures (5.26) and (5.27). It is observed 
that the displacement along the X-direction is significantly lesser for the frictional case, 
more so at the bottom node. Further, in the frictional case, the ^-displacement plots 
of centre and bottom nodes are very much different, while they are almost same in the 
frictionless case. This is due to the rotation of the sphere induced by the frictional forces, 
which reduce the displacement along X-direction, but increase the displacement along Z- 
direction at the bottom node. This rotation of the sphere in the frictional case is clearly 
coon in fiomrp<5 ( F, 99 - Fi 9.5V Alfio. it’is observpfi that the disnlacement alone negative 
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Figure 5.18: Impact of sphere on plate (Prictionless): Configuration at t 0.0 

The velocity given to the sphere is (10,0,-10) m/s, in the coordinate system shown in 
figure (5.15) and the coefiicient of friction at the contact interface is 0.3, for the case of 
analysis with friction. The time step used in this analysis is 0.0003 seconds. The X and 
Z displacements of the bottom node and the centre of the sphere for both the frictionless 
and the frictional cases are shown in figures (5.29) and (5.30). The variation of maximum 
equivalent stress in both the sphere and plate for both cases is shown in figure (5.31)- 
Convergence difficulties are observed in the frictional case, and the analysis could be done 
till 0.0123 seconds, while the contact between the bodies are lost at about 0.02 seconds. 
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5.25: Impact of sphere on plate (Frictional): Configuration at t = 0.06 



Figure 5.26: Displacement of the sphere along X-direction 
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Figure 5.27; Displacement of the sphere along Z-direction 
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Figure 5.28: Maximum equivalent stress in sphere and plate 
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Figure 5.29; Displacement of the sphere along X-direction (polycarbonate material) 



Figure 5.30: Displacement of the sphere along Z-direction (polycarbonate material) 
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Figure 5.31: Maximum equivalent stress in sphere and plate (polycarbonate material) 




Chapter 6 


Conclusions and Scope for Future 
Work 


6.1 Conclusions 

A finite element code is developed for solving 3-D, dynamic, elasto-plastic contact prob- 
lem including the effects of friction at the contact interface. A node-to-segment contact 
formulation is used for developing the contact stiffness matrix. The Lagrange multiplier 
method is used for enforcing the contact constraints. Frictional effects are included using 
the classical Coulomb’s friction law. The updated Lagrangian formulation along with 
the modified Newton-Raphson iterative scheme is used for carrying out the elasto-plastic 
analysis. Divergence handling procedures, viz., under-relaxation and line search methods 
are incorporated to improve the convergence characteristics of Newton-Raphson itera- 
tions. Green-Lagrange strain tensor and objective Jaumann stress measure are used in 
the formulation of the governing equations. Elasto-plastic behaviour is modelled by an 
associated flow rule based on von Mises yield criterion and power law type isotropic hard- 
ening. Euler forward integration is used for the integration of the constitutive equation. 
Newmark’s method, which is the best finite difference method in terms of stability and 
accuracy, is used for the time integration. An unloading scheme is incorporated to inlude 
the effect of local or global unloading. 

The code is first validated by solving some standard problems, for which, the results 
are available in the published literature. Then one additional problem, the elasto-plastic 
impact of sphere on axi elasto-plastic plate is solved to demonstrate the applicability of 



the program. Therefore, this program can be used for analysing the practical problems. 


6.2 Scope for Future Work 

The desired future extensions of this work are the following: 

1. Currently the program uses Green-Lagrange strain measure and Jaumann stress 
measure. Other strain and stress measures can be incorporated and comparative 
study of the results can be performed. 

2. Other types of finite elements can be incorporated. 

3. Material damping effects can be included. 

4. Option for using penalty function method for developing the contact stiffness matrix 
can be incorporated. 

5. Other material models like kinematic and anisotropic hardening and hardening laws 
other than power law can be implemented. 

6. Non-classical friction laws can be used for modelling the effect of friction. 

7. Finally, it is highly desirable to have versatile pre and post processors, so that the 
program can be more user-friendly. 
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